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ABSTRACT 

Detailed theoretical studies of the high-redshift universe, and especially reionization, are generally 
forced to rely on time-consuming N-body codes and/or approximate radiative transfer algorithms. We 
present a method to construct semi- numerical "simulations" , which can efficiently generate realizations 
of halo distributions and ionization maps at high redshifts. Our procedure combines an excursion- 
set approach with first-order Lagrangian perturbation theory and operates directly on the linear 
density and velocity fields. As such, the achievable dynamic range with our algorithm surpasses 
the current practical limit of N-body codes by orders of magnitude. This is particularly significant in 
studies of reionization, where the dynamic range is the principal limiting factor because ionized regions 
reach scales of tens of comoving Mpc. We test our halo-finding and ionization-mapping algorithms 
separately against N-body simulations with radiative transfer and obtain excellent agreement. We 
compute the size distributions of ionized and neutral regions in our maps. We find even larger ionized 
bubbles than do purely analytic models at the same volume- weighted mean hydrogen neutral fraction, 
xhi, especially early in reionization. We also generate maps and power spectra of 21-cm brightness 
temperature fluctuations, which for the first time include corrections due to gas bulk velocities. We 
find that velocities widen the tails of the temperature distributions and increase small-scale power, 
though these effects quickly diminish as reionization progresses. We also include some preliminary 
results from a simulation run with the largest dynamic range to date: a 250 Mpc box that resolves 
halos with masses M > 2.2 x lO^M©. We show that accurately modeling the late stages of reionization, 
^Hi ^0-5, requires such large scales. The speed and dynamic range provided by our semi- numerical 
approach will be extremely useful in the modeling of early structure formation and reionization. 
Subject headings: cosmology: theory - early Universe - galaxies: formation - high-redshift - evolution 



1. INTRODUCTION 

Accurately modeling the formation of bound struc- 
tures is invaluable for understanding any process in 
the early universe. Reionization, the epoch when 
radiation from early generations of astrophysical objects 
managed to ionize the intergalactic medium (IGM), is 
particularly sensitive to the distribution of collapsed 
structure. Current observations p aint a complex pic- 
ture of the reionizatio n epoch (iMesinger fc Haima: 
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2004; 



20041: IWvithe fc Loebj l2004l: iFan et al.. 
Mesinger HaimanI l2006l : iMalhotra fc R.hoadsl 
Furlanetto et al." "2006c; "Malhotra & Rhoads" 2006; 
Page et al. 2006; Kashikawa et al. 2006; Totani et al. 
2006). The next generation of instruments (James 
Webb Space Telescope; 21-cm instruments such as the 
Low Frequency Array and the Mileura Widefield Array 
Low- Frequency Demonstrator; CMB polarization mea- 
surements with Planck^ etc.), could potentially shed light 
on this poorly understood milestone. Unfortunately, 
we still do not have accurate models of reionization 
with which to interpret these upcoming {and current) 
observations. 

The main difficulty lies in the enormous dynamic range 
required. Ionized regions are expected to reach charac- 
teristi c sizes of tens of comov ing Mpc ( Furlanetto et al. 
[2004^ iFurlanetto fcOhl[2005h . which is over seven or- 
ders of magnitude in mass larger than the pertinent cool- 
ing mass, correspo nding to g as wit h a temperature of 
T ^ 10^ K (e.g. [Ife tathimJ ^992"; "ThouLAWeinberg 
[T996[ lGnedidl2000l3 : [Shapi ro et al. 1991). The required 
dynamic range is even larger if smaller "minihalos" be- 



low this cooling threshold are important during reion- 
ization. Because of the steep mass dependence of halo 
abundances, halos with masses close to the cooling mass 
could dominate the photon budget. Hence modeling 
reionization requires simulation box sizes of hundreds 
of megaparsecs on a side, with extremely high resolu- 
tion. Attempts to overcome these obstacles have gener- 
ally followed the same fundamental and well-trod path 
(e.g. Gnedin 2000a; Razoumov et al. 2002; Ciardi e t al. 
[2003: Sokasian et al. 2003; Iliev et al. 2006b; Zahn etHI 
120071 : LTrac & Cen 2006): (1) N-body codes are run to 
generate halo distributions; (2) a simple prescription is 
used to relate the halo mass to an ionizing efficiency; (3) 
approximate methods (generally so-called ray-tracing al- 
gorithms) are used to model radiative transfer (RT) on 
large scales. 

Even with modest halo resolution 

(jSpringel fc HernquistI l2003l ) of tens of dark matter 
particles per halo, such schemes are computationally 
limited to box sizes of tens of mega parcecs, if they 
wish to resolve the likely cooling mass. McOuinn e t al.l 
(2006a) extended the mass resolution of their sim- 
ulations by using a merger tree scheme to populate 
sub-grid scales with unresolved halos in a stochastic 
manner. Such hybrid schemes are useful for extending 
the dynamic range, but merger trees require a number 
of c orrections to achieve consistent mass funct ions (see, 
e.g. ISheth PitmanI [19971 : iBenson et al.l l2005L and Fig. 
1 in lMcOuinn et alll2006a ) and to track individual halos 
with redshift. Moreover, although they are perfectly 
adequate for many purposes (including studying the 
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large-scale features of reionization), they prevent one 
from taking full advantage of the simulation. 

Aside from dynamic range, the other main limiting fac- 
tor in all of the above numerical approaches is speed. 
Even if the relevant scales can be resolved with N-body 
codes, such as may be the case in the early phase of reion- 
ization or with hybrid stochastic schemes. The codes 
themselves generally take days to run on large super- 
computing clusters, with the approximate RT algorithms 
consuming a few additional days. The computational 
cost of each simulation makes it difficult to explore the 
full range of parameter space for reionization, which is 
particularly large because we know so little about high- 
redshift galaxies. 

The computational cost becomes truly prohibitive if 
hydrodynamics is included: the largest such simulation 
of rei onization performed to date spanned only 10h~^ 
Mpc (|Sokasian et al.|[2QQ3ll ). Including self-consistent de- 
scriptions of galaxy formation - even at the approximate 
level currently implement ed in low er-redshift cosmolog- 
ical simulations (e.g., Sp ringel fc H ernquist 2003) - re- 
quires hydrodynamics, so N-body simulations of reion- 
ization are limited to semi-analytic prescriptions for star 
formation, feedback, etc. It is therefore worthwhile to 
explore even simpler schemes. 

The purpose of this paper is to introduce approximate 
but efficient methods for generating halo distributions at 
high redshifts as well as for generating the associated ion- 
ization maps. W e apply an excursion -set approach (e.g. 
[Bond et alJliMil; iLacev Colel ll993[ ) to the filtering of 
a realization of the linear density field and then adjust 
halo locations with first-order perturbation theory. We 
can thus generate halo distributions at any given red- 
shift, without explicitly including information from any 
higher redshifts. This scheme is an updated form of 
the ''peak -patch" formalism developed and validated by 
iBond fc MversI (|l996al lbl). although it was conceived and 
implemented completely independently. We then apply 
a similar technique to obtain the ionization field from the 
ha lo field. Thi s part is similar to the schemes described 
in IZahn et al.l ([2005, 2007), except applied to our effi- 
ciently built halo distributions. As such, our methods 
allow us to make general predictions about non-linear 
processes, such as structure formation and reionization, 
without making use of time-guzzling cosmological sim- 
ulations. The speed of our approach also allows us to 
explore a larger dynamic range than is possible with cur- 
rent cosmological simulations while preserving detailed 
spatial information (at least in a statistical sense), un- 
like purely analytic models. 

This paper is organized as follows. In §[21 we introduce 
and test the components of our halo finding algorithm. 
In § [3l we introduce and test our HII bubble finding 
algorithm. In § HI we use our semi-numerical scheme 
to generate maps and power spectra of expected 21-cm 
brightness temperature fluctuations throughout reioniza- 
tion. In § O we summarize our key findings and present 
our conclusions. 

Unless stated otherwise, we quote all quantities in co- 
moving units. We adopt the background cosmological 
parameters (I^a, flu, ^6, ^, crs, ^o) = (0.76, 0.24, 
0.0407, 1, 0.76, 72 km s"^ Mpc"^), consistent wit h the 
three -year results of the WMAP satellite (S per gel et al.l 
[2006h . 



2. SEMI-NUMERICAL SIMULATIONS OF HALO 
PROPERTIES 

In brief, our algorithm generates a linear density field 
and identifies halos within it. Because only linear evo- 
lution is required, the algorithm is fast and flexible. We 
generate 3D Monte-Carlo realizations of the linear den- 
sity field on a box with sides of length L = 100 Mpc and 
N = 1200^ grid cells. As such, we are able to take advan- 
tage of many pre-existing tools operating on the linear 
density field alone. Our method consists of the following 
principal steps: 

1. creating the linear density and velocity fields 

2. filtering halos from the linear density field using 
the excursion-set formalism 

3. adjusting halo locations using their linear-order 
displacements 

Step (1) only needs to be done once for each realization, 
since it is independent of redshift. As mentioned above, 
steps (2) and (3) need only be performed on redshifts of 
interest, i.e. since our output at redshift z is independent 
of any outputs at higher redshifts, there is no need for 
our code to "run down" to as is the case for N-body 
codes. 

Our algorithm is an updated and simplified version o f 
the "peak-patch" algorithm of iBond fc MyersI (|l996af ): 
we refer the interested reader there for more detailed 
explanations of some steps. A simpler version has also 
been used by Scannapieco et al. (2002) to study metal 
enrichment at high redshifts. 

We perform our semi-numerical simulations on a sin- 
gle desktop Mac Pro with two dual-core 3.00 GHz Quad 
Xeon processors and 16 GB of RAM. For N = 1200^, step 
(1) takes ~ 1 hour. For a given redshift in our range of 
interest, specifically for z = 8.75, steps (2) and (3) take 
~ 2.5 hours. To achieve comparable halo mass resolution 
(including halos with M > lO^M©) with a minimum of 
~ 500 particles per halo ([Springel fc Hernqui st "20031. N- 
body codes would require a prohibitively large number 
of particles, N ~ 10^^! Below we describe in detail the 
components of our model. 

2.1. The Linear Density Field 

Our linear density field is generated in much the same 
way as it is for N-body codes. We briefly outline the 
procedure here. 

The density field of the universe, 5(x) = p(x)/p — 1, 
in the linear regime^ is well-represented as a Gaussian 
random field, whose statistical properties are fully de- 
fined by its power spectrum, cr^(k) = (|5(k)p). Here, 
J(k) is the Fourier transform of ^(x), and the standard 
assumption of isotropy implies cr^(k) = cr^(/c) while ho- 
mogeneity implies that there are no density fluctuations 
with wavelengths larger than the box size L = yi/^. We 
use the followin g;, standard (e.g. IBagla fc PadmanabhanI 
I1997I : ISirkd[2QQ5.) Fourier transform conventions: 

S{k) = ^^S{^)e-''^-- , (1) 

^ In linear theory, density perturbations evolve in redshift as S(z) 
= 6(0) D(z), where D{z) is the linear growth factor normalized 
so that D(0) = 1 (e.g., [Tiddle et al. 1996). Unless the redshift 
dependence is noted explicitly, from this point forward we will work 
with quantities linearly-extrapolated to z = 0. 
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(2) 



with the inverse transform being 

The discrete simulation box only permits a finite set 
of wavenumbers: k = A/c(z, j, /c), where Ak = 27t/L and 
i, j, /c are integers in the range {-\^N /2, /2]. For each 
independent wavenumber,^ we assign 



(5(k) 



(3) 



where ak and 6k are drawn from a zero-mean Gaussian 
distribution with unit varian ce. We use the power spec- 
trum from Eisenstein fc Hu| |T999). Then the real-space 
density field, (^(x), is obtained by performing an inverse 
Fourier transform on ^(k). 

2.2. The Linear Velocity Field 

We construct a linear velocity field corresponding to 
our linear density fi eld usin^ the standard Zel'Dovich 
approximation ( c.f. IZerPovichl ll97Ql : lEfstathiou et al.l 
[iMl;[Sir3[2005|): 



(4) 
(5) 



v = xi V^(x) , 

(^(x) = -V-[(l + (^(x))^(x)] 

^-V-V^(x) , (6) 

where x and xi denote initial (Lagrangian) and updated 
(Eulerian) coordinates, respectively, V^(x) is the displace- 
ment vector, and the last equation follows from the conti- 
nuity criterion, with the final approximation using linear- 
ity, ^(x) <C 1. We note again that all units are comoving, 
unless stated otherwise. From the above, one can relate 
the velocity mode in our simulation at redshift z to the 
linear density field: 



(7) 



where for computational convenience differentiation is 
performed in /c-space. 

Another convenient property of this first-order 
Zel'Dovich approximation is that the velocity field can 
be decomposed into purely spatial, Va;(x), and purely 
temporal, Vz{z), components: 

v(x, = Vz{z)vx{-x.) , (8) 

where Vz{z) = D{z) and Vx{x.) is the inverse Fourier 
transform of ik^(k)//c^. This is computationally conve- 
nient, as we only need to compute the Vx{x.) field once 
in order to be able to scale it for all redshift s, and it 
also allows us to write a simple, exact expression for the 
integrated linear displacement field, ^. When eq. ([8]) is 
integrated from some large initial zq [D{zo) <C D{z)]^ the 
total displacement is just 



*(x,^) = [D(^)-Z)(zo)]v,(x) 
«i?(^K(x) 



(9) 



^ Since 6{x.) is real-valued, only half of the A^-modes defined 
above are independent. The other half are determined by the usual 
Hermitian constraints for r eal- valued functions (see for example 
IHocknev Eastwood! [19881 : IBagla PadmanabhanlfT997h . 



We make use of this displacement field to adjust the halo 
locations obtained by our filtering procedure (see § 12. 4p . 
as well as to adjust the linear density field for our 21-cm 
temperature maps (see §|4]). 

In principle, one could obtain non-linear velocities by 
mapping the linear overdensity to a corresponding non- 
linear Qverdensity obtaine d from a spherical collapse 
model (|Mo fc White! Il99"6l ). and then taking the time 
derivative of the non-linear overdensity. However, due 
to the large spread in the dynamical times of the non- 
linear density field, accurately capturing the time evolu- 
tion is non-trivial. Furthermore, although the non-linear 
density field implicitly captures the velocities of collaps- 
ing gas, mapping each pixel's linear density to its non- 
linear counterpart independently of other nearby pixels 
does not properly preserve correlations on larger scales. 
Hence, we choose to use the linear density field directly 
in estimating velocities. For the purposes of studying the 
ionization field, we are further justified in this procedure 
because our final ionization maps are smoothed on large 
scales, on which most pixels are still in the linear regime 
at the high redshifts of interests. It is possible to include 
higher-order contributions to the ZePDovich approxini a- 
tion where necessary (e.g.. !Scoccimarro fc Shet'hl!2QQ2f ). 

2.3. Halo Filtering 

In standard Press- Schechter theory (PS; see e.g.. 
Press fc Schechted!l974l : !Bond et al.![l99l! : !Lacev fc Cold 
1993! ). the halo mass function can be written as 



dn{> M, z) 
dM 



2 p Sc{z) da{M) 



■ exp 



2(j2(M) 
(10) 

where n{> M, z) is the mean number density of halos 
with total mass greater than M, p = r^MPcrit is the 
mean background matter density, Sc{z) ~ 1.68/ D{z) is 
the scale-free critical over-density evaluated in the case 
of spherically symmetric collapse ([Peebles! !l98Qf ). and 
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Jfdk 



cj\k)W\k,M) , (11) 



is the squared r.m.s. fiuctuation in the mass enclosed 
within a region described by the filter function, W{k^ M), 
normalized to integrate to unity. 

Although the PS mass function in eq. (p!Q|) is in fair 
agreement with simulations, especially for halos near the 
characteristic mass, at low redshifts it underestimates 
the number of high-mass halos and overestimates the 
number of low-mass halos when compared wit h large nu- 
merical simulations (e.g. ! Jenkins et al.!!2QQll ). A mod- 
ified expression shown to fit low-redshift simulation re- 
sul ts more accurate l y (to within ~ 10%) was obtained 
bv lSheth fc Tormen! ([T999h : 



dn{> M,z) 



p d[\n (j{M)] 
M dM 



TT 



1 

0^ 



ve 



(12) 

where z> = \ /a6riz) /a(M \ and a, p, and A are fitting 
parameters. !Sheth et al.! (pOOl) derive this form of the 
mass function by including shear and ellipticity in model- 
ing non-linear collapse, effectively changing the scale- free 
critical over-density Sc{z)^ into a function of filter scale. 



4 



5c{M, z) = VaSc{z) 



aHM)Y 



(13) 



Here b and c are additional fitting parameters (a is 
the same as in eq. [T2]) . For the consta nts above, w e 
adopt the recent values obtained by Jenkins ^et alJ (|2QQll ). 
who studied a large range in redshift and mass: a = 
0.73, A = 0.353, p = 0.175, b = 0.34, c = 0.81. We 
note, however, that the situation at high redshift s is 
less clear: studies disa gree o n the relative accur acy o f 
the Press & Schechte il (Il97l and I Jenkins et all (120011) 
form s ([Reed et alj 120031 : llliev et al.l l2006bl : IZahn et all 
|2007[ ). Our algorithm can be trivially modified to ac- 
commodate other choices for the mass function; fortu- 
nately, for the purposes of the ionization maps (see ^Sj), 
the choice of mass function makes very little difference 
because all have a similar de pendence on the local den- 
sity (Furl anetto et al.lL2006al ). 

The mass functions in equations (p!Q|) and (p!2|) can be 
obtained by the standard excursion set random walk pro- 
cedure. The approach is to smooth the density field 
around a point, x, on successively smaller scales start- 
ing with M ^ oo [where cr^(M) 0] and to identify 
the point as belonging to the halo with the largest M 
such that (^(x, M) > 6c{M, z). If M) is chosen to 

have a sharp cut-off, this procedure amounts to a random 
walk of J(x, M) along the mass axis, since the change in 
J(x, M) as the scale is shrunk is independent of (5(x, M) 
for a top-hat filter in /c-space (see eq. [TT]) . 

We perform this procedure on our realization of the 
linear density field by filtering the field using a real- 
space top-hat filter^, starting on scales comparable to 
the box size and going down to grid cell scales, in log- 
arithmic steps of width AM/M = 1.2."^ At each filter 
scale, we use the scale-dependent barrier in eq. ([13]) 
to mark a collapsed halo if (5(x, M) > 5c{M^z). Filter 
scales large enough that collapsed structure is extremely 
unlikely, 5c{M^z) > 7<j(M), are skipped (Mesinger et all 
l2005i) . Since this procedure treats each cell as the center 
of a spherical filter, neighboring pixels are not properly 
placed in the same halo. Because of this, we discount 
halos which overlap with previously marked halos. 

As mentioned above, this algorithm is similar 
to the ^^peak - patch" approach first introduced by 
[Bond fc Mversl (|l996af ). The primary differences are: 

^ There is a slight swindle in the current application of this for- 
malism. The filter function is assumed to be a top-hat in k-space 
in order to facilitate the analytic random walk approach described 
above. However, when the power spectrum is normalized to ob- 
servations [i.e. (j(R = 8/i~-'^Mpc] = as), the filter that is used to 
define the mass M corresponding to is a top-hat filter in real 
space. Nevertheless, it has been shown that the mass function is 
not very sensitive to this filter choice (Bond et al. 1991). 

^ We note that Mesinge r et al .| (|2005l) required a much smaller 
step size at these redshifts, AM/M ~ 0.1, in order to produce ac- 
curate mass functions using ID Monte-Carlo random walks. How- 
ever, here we find that we can reproduce accurate mass functions 
with a larger step size, since in our 3D realization of the density 
field, "overstepping" 6c{M, z) due to a large filter step size can be 
compensated with a small offset in the filter center, i.e. by cen- 
tering the filter in a neighboring cell. This is the case since over- 
stepping (5c (M, z) means that some dense matter between the two 
filter scales was "missed" . In a ID Monte-Carlo random walk this 
matter is unrecoverable; however, in a 3D realization of the density 
field, the missed matter will be picked up by a filter centered on a 
neighboring cell. 



(1) we use the I Jenkins et al.l (|2QQll ) barrier to identify 
halos (rather than calculating the strain tensor to ac- 
count for ellipsoidal collapse), (2) we do not separately 
identify peaks in the density field (this step is not re- 
quired given modern computing power), and (3) we use 
the ^^full exclus i on" cri terion for preventing halo overlap. 
iBond fc Mversl (|l996bl ) found that a "binary exclusion" 
method in which pairs of overlapping halos are compared 
and eliminated was somewhat more accurate. However, 
at the high redshifts of interest to us, halo overlap is rare, 
and we are primarily interested in the large-scale prop- 
erties of the halo field, which are relatively insensitive to 
the details of the overlap criterion. 

We also note that our halo finder is similar 
in spirit to the PT Halos algorithm introduced by 
IScoccimarro fc ShethI ([2002) to generate mock galaxy 
surveys at low redshifts. There are two key differences. 
First, at present we use only first-order perturbation the- 
ory to displace the particles.^ This limits us to higher 
redshifts, where velocities are smaller. However, our al- 
gorithm does not require particles in order to resolve ha- 
los and hence can accommodate a considerably larger 
dynamic range than PTHalos. 

Mass functions resulting from this procedure are shown 
as points in Figure [H with error bars indicating l-a Pois- 
son uncertainties and bin widths spanning our mass filter 
steps. Dotted red curves denote PS mass functions gen- 
erated by eq. (p!Q|) : short-dashed blue curves denote ex- 
tended PS conditional mass functions generated by eq. 
([To]) but also taking into account the absence of den- 
sity modes longer than the box size; long-dashed green 
curves denote mass functions generated using the Sheth- 
Tormen correction in eq. ([T2|) . The upper (lower) set 
of curves and points correspond to redshifts of z=6.5 
(z=10). The dotted and short-dashed curves overlap at 
these redshifts due to our large box size (L = 100 Mpc), 
so we are immune to th e finite box effects pointed out by 
iBarkana fc Loebl (|2004[ ). 

Fig. [U shows that we obtain accurate mass functions 
for M > lO^M©. Our procedure seems to underpre- 
dict the abundance of halos with masses approaching 
the cell size, Mceii ~ lO^M©. However, as the Jeans 
mass corresponding to a gas temperature of ~ 10^ K 
is Mj{z ~ 8) ~ IO^Mq, in subsequent calculations, we 
only use halos with masses greater than Mmin = lO^M©. 
Using this Mmin, we match the collapse fraction obtained 
by integrating eq. ([T2|) to better than ~ 10%. 

This mass cutoff corresponds to the minimum tem- 
perature required for efficient atomic hydrogen cooling 
and would be the pertinent mass scale if: (1) the H2 
cooling channel is suppressed, e.g. due to a perva- 
sive Lyman- Werner (LW) background, and if (2) photo- 
ionization feedback is ineffective at suppressing gas cool- 
ing and collapse onto higher mass halos. While feed- 
back at high redshifts remains poorly-constrained, both 
of these assumptions seem reasonable during the mid- 
dle stages of reionization on which we focus. A dis- 
sociating LW background is likely to have established 

^ We note that a similar scheme to ours has been independently 
created by O. Zahn (private communication). This scheme uses 
a simple Press-Schechter barrier but adjusts halo locations follow- 
ing second-order Lagrangian perturbation theory. However, he has 
found that the second-order corrections make very little difference 
to the map. 
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Fig. 1. — Mass fu nctio ns generated from our halo filtering pro- 
cedure discussed in ^2.31 are shown as points. Dotted red curves 
denote PS mass functions generated by eq. (I10|) : short-dashed 
blue curves denote extended PS conditional mass functions gener- 
ated by eq. (|10|) but also taking into account the absence of density 
modes longer than the box size; long-dashed green curves denote 
mass functions generated using the Sheth-Tormen correction in eq. 
(|12|) . The upper (lower) set of curves and points correspond to 
redshifts of z=6.5 (z=10). 



itself well before th e universe is significantly ionized 
(|Haiman et alJ [19971). Model-dependent empirical evi- 
dence supporting the suppression of star formation in 
smaller mass halos, M < Mmin, can also be gleaned 
from WMAP data (|Haiman fc BrvanI I2QQ6D . Further- 
more, although early work suggested that an ionizing 
background could partially suppress star formation in 
halos with virial temperatures of Tvir < 3.6 x 10^ K 
(M < 2 xlQ^M0) dThouLfc Weinberg 1996), more recent 
studies (|Kitavama fc Ikeuch i 2000 : D iikstra et al. 2004) 
find that at high redshifts {z > 3), self-shielding and the 
increased cooling efficiency could be strong countering 
effects for halos with virial temperatures Tvir > 10^ K. 
We postpone a more detailed analysis of the reionization 
footprint left by photo-ionization feedback to a future 
work. 

2.4. Adjusting Halo Locations 

Once the halo field is obtained, we use the displace- 
ment field obtained through eq. (|9|) to adjust the halo 
locations at each redshift. This corrects for the enhanced 
halo bias in Eulerian space with respect to our filtering, 
which is done in Lagrangian space (i.e. using the initial 
locations at large z). For computational convenience, we 
smooth the 1200^ velocity field onto a coarser-grained 
200^ grid before adjusting halo locations. The choice of 
resolution, where each cell is (100 Mpc)/200 = 0.5 Mpc 
on a side, is somewhat arbitrary here, and we have veri- 
fied that our halo and 21-cm power spectra are unaffected 



Fig. 2. — Halo power spectra at z = 8 .7, with L = 20 / i~^Mp c 
and cosmological parameters taken from IMcQuinn et al.l (|20Q6al V 
The solid red curve is the halo power spectrum from an N-body 
simulation obtained from McQuinn et al.l (|2006al) (c.f. the bottom 
panel of their Fig. 2). The short-dashed green and the long-dashed 
violet curves are obtained from our filtering procedure with and 
without the halo location adjustments, respectively. 

by this choice. We also note that in linear theory, the 
mean velocity dispersion inside a (0.5 Mpc)^ sphere with 
mean density at 2: = 10 is a factor of ^10 lower than the 
r.m.s. bulk velocity of such regions, so smoothing over 
smaller scale velocities appears reasonable. Furthermore, 
we keep in mind that our "endproducts" in this work are 
ionization and 21-cm temperature ffuctuation maps, for 
which such "low-resolution" is more than adequate (com- 
pare, e.g., to N-body simulations of reionization, which 
typically have similar cell sizes for the radiative transfer 
component). 

In Figure O we plot the halo power spectrum, de- 
fined as Ahh(^,^) = A:^/(27r^V) (|5hh(k, 2;)^)/^, where 
4h(x, ^) = Mcoii(x,z)/(Mcoii(^)) - 1 is the collapsed 
mass field. ^ The solid red curve is the halo power spec- 
trum from a 2 Mpc N-body sim ulation at z = 8.7 
obtained from iMcQuinn et_alj (j2006al ) (c.f. the bottom 
panel of their Figure 2). The short-dashed green and the 
long-dashed violet curves are obtained from our filtering 
procedure (matching the assumed cosmology) with and 
without the halo location adjustments, respectively. We 
note that ignoring the cumulative motions of halos re- 
sults in an underestimate of the power of long- wavelength 
modes of the halo field by a factor of ~ 2 in this case. 
The average Eulerian bias of these halos is ~ 2, about 
half of which comes from the correction from Lagrangian 
to Eulerian coordinates. 

After the halo locations are adjusted according to lin- 
ear theory, our halo power spectrum agrees almost per- 

^ We use the collapsed mass field, rather than the individual 
galaxies, because we calculate the power from the smoothed cells. 
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fectly with the simulation. By design our procedure in- 
cludes Poisson fluctuations in the halo number counts, 
which dominate the power spectrum at > 5 /i/Mpc 
and ar e lost in purely analytic estimates (jMcQuinn et alJ 
l2QQ6a| ). We also note that both the halo mass func- 
tions and power spectra are statistical tests and hence 
the agreement shown here does not imply that our halo 
field has a one-to-one mapping with an N-body halo 
field so urced by i dentical initial conditions. Indeed, 
[Gelb fc Bertschingerl (|l994f ) showed that those particles 
located nearest initial linear density peaks are not nec- 
essarily incorporated into massive galaxies. The "peak 
particle" algorithm is less robust than our smoothing 
technique, but we still do not expect to recover halo 
masses or locations precisely. We plan on doing a "one- 
on-one" comparison between halo fields obtained from 
our halo finder to those obtained from N-body codes 
in a future work. However, it is certainly encouraging 
that t he very similar " peak-patch" group finding formal- 
ism of lBond (p^96a) did very well when com- 
pared "one-on-qne" to N-body codes at large mass scales 
(|Bond M^ll996br ). 

In Figure [3] we show slices through the halo field from 
our simulation box at z = 8.25, generated by the above 
procedure, again with {right panel) and without {left 
panel) the halo location adjustments. In the figure, the 
halo field is mapped to a lower resolution 400^ grid for 
viewing purposes. Each slice is 100 Mpc on a side and 
0.25 Mpc deep. Collapsed halos are shown in blue. Vi- 
sually, it is obvious that peculiar motions increase halo 
clustering. 

3. GENERATING THE IONIZATION FIELD 

Once the halo field is generated as described above, 
we can perform a similar filtering procedure (also us- 
ing the excursion-set formalism) to obtain th e ionization 
field (similar methods have been discussed bv lZahn et al.l 
[20051 . 12007i) . The time required for this final step is a 
function of xhi, with large xhi requiring less time than 
small xui- Specifically, at xm ~ 0.5 this step takes ~ 15 
minutes to generate a 200^ ionization box on our work- 
station. 

There are two main differences between the halo fil- 
tering and the HII bubble filtering procedures: (1) HII 
bubbles are allowed to overlap, and (2) the excursion 
set barrier (the criterion for ionization) becomes, as per 
iFurlanetto et all (|2004af ): 

/coii(xi,M,^)>r' , (14) 

where ( is some efficiency parameter and /coii(xi, M, z) 
is the fraction of mass residing in collapsed halos inside 
a sphere of mass M = 4/37ri?^p[l + {Sn\{'^i^ z))r\^ with 
mean physical overdensity ((5ni(xi, z))r^ centered on Eu- 
lerian coordinate xi, at redshift z. 

Equation (p!4|) is only an approximate model and makes 
several simplifying assumptions about reionization. In 
particular, it assumes a constant ionizing efficiency per 
halo and ignores spatially-dependent recombinations and 
radiative feedback effects . It can easily b e modified to 
include th ese effects f e.g.. IFurlanetto et al.p 2004b. 2006a; 
IFurlanett o fc: Ohl ii2005). and we plan to do so in future 
work. Here we present the simplest case in order to best 
match current RT numerical simulations. 



This prescription models the ionization field as a two- 
phase medium, containing fully-ionized regions (which 
we refer to as HII bubbles) and fully- neutral regions. 
This is obviously much less information than can be 
gleaned from a full RT simulation, which precisely tracks 
the ionized fraction. However, HII bubbles are typi- 
cally highly-ionized during reionization, and for many 
purposes (such as for 21 cm maps) this two-phase ap- 
proximation is perfectly adequate. 

In order to "find" the HII bubbles at each redshift we 
smooth the halo field onto a 200^ grid. Then we filter 
the halo field using a real-space top-hat filter, starting on 
scales comparable to the box size and decreasing to grid 
cell scales in logarithmic steps of width AM/M = 0.33. 
At each filter scale, we use the criterion in eq. (p!4|) 
to check whether the region is ionized. If so, we flag 
all pixels inside that region as ionized. We do this for 
all pixels and scales, regardless of whether the resulting 
bubble would overlap with other bubbles. Note, there- 
fore, that the nominal ionizing efficiency ( that we use 
as an input parameter does not equal (1 — xhi)//co11- 
They typically differ by < 30%, with C/coii < 1 - ^Hi 
early in reionization and C/coii > 1 — ^hi late in reioniza- 
tion). Unfortunately, we thus cannot use our algorithm 
to self-consistently predict the time evolution of the ion- 
ized fraction (rather, that must be prescribed from some 
other model). Of course, the same is true for N-body 
simulations, because the evolution of the ionized fraction 
depends on the evolving ionization efficiency of galaxies 
and cannot be self-consistently included in any present- 
day simulation. 

In order to obtain the density field used in eq. (p!4|) . 
^^ni(xi,z), we use the Zel'Dovich approximation on our 
linear density field, ^(x), in much the same manner as 
we did to adjust our halo field in § 12. 4[ Starting at some 
arbitrarily large initial redshift (we use zq = 50), we 
discretize our high-resolution 1200^ field into "particles" 
whose mass equals that in each grid cell. We then use 
the displacement field (eq. [9|) to move the particles to 
new locations at each redshift. This resulting mass field 
is then smoothed onto our lower resolution 200^ box to 
obtain (^ni(xi,z). We then recalculate the velocity field 
(§ 12.21) using t he ne w densities. 

Za hn et al.l (|2007l ) showed that a very similar HII bub- 
ble filtering procedure performed on an N-body halo field 
was able to reproduce the ionization topology obtained 
through a ray-tracing RT algorithm fairly well. Their 
algorithm differs from ours in two ways. First, they used 
a slightly different barrier definition; however, this dif- 
ference has only a small impact on the ionization topol- 
o^yJ More importantly, for each filter scale at each pixel, 
IZahn et"al] ([200 7) flag only the center pixel as ionized if 
the barrier is crossed, whereas we flag the entire filtered 
sphere. 

In order to test our bubble filtering algorithm, we ex- 
ecute it on the same N-body halo field at z = 6.89 as 
was used t o gene rate the bottom panels of Fig. 3 in 
IZahn et al.1 ^U^. We compare analogous ionization 
maps created using various algorithms in Figure [H All 

^ Specifically, in order to match the physics of their simulations 
better, they required / fit /coll > C""*^- However, the density mod- 
ulation ends up nearly identical to our model, so the topology is 
almost unchanged. 
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Fig. 3. — Slices through the halo field from our simulation box at z — 8.25. The halo field is generated on a 1200^ grid and then mapped 
to a 400^ grid for viewing purposes. Each slice is 100 Mpc on a side and 0.25 Mpc deep. Collapsed halos are shown in blue. The left panel 
shows the halo field directly filtered in Lagrangian space; the right panel maps the field to Eulerian space according to linear theory (see 
^ 12.41 and eq. [9]). The right panel corresponds to the bottom-left (xhi = 0.53) ionization field in Figure[5] 



slices are 93.7 Mpc on a side and 0.37 Mpc deep, with ^ 
adjusted so that the mean neutral fraction in the box is 
^Hi = 0.49. Ionized regions are shown as whi te. The left- 
most and right-most panels are taken from IZahn et al] 
(j2007D . The left-most panel was created by perform- 
ing their bubble filtering procedure directly on the lin- 
ear density field (without explicitly identifying halos). 
The second panel was created by performing their bub- 
ble filtering procedure on their N-body halo field, but 
with the slightly different barrier definition in eq. 
The third panel was created by performing our bubble 
filtering procedure on the same N-body halo field, but ig- 
noring density fiuctuations outside of halos (i.e. setting 
((5ni(xi, = 0), which we have verified give nearly 
identical bubble maps as our full procedure (so long as 
xhi is fixed). The right-most p anel w as created using 
an approximate RT algor ithm (|Abel fc Wandelt 2002; 
iSokasian et 111120011 . [20031 ) on the same halo field. 

It is immediately obvious from Fig. 2] that all of the 
approximate maps {first three panels) reproduce the RT 
map {right-most panel) fairly well. Even the HII bub- 
ble filtering performed directly on the linear density field 
{left-most panel) performs well, which is encouraging, as 
that is the starting point for our semi-numerical proce- 
dure and we only improve on this scheme. 

Figure m shows that our HII bubble filtering algorithm 
is an excellent ap proximation to RT . The similar algo- 
rithm proposed bv lZahn et al.l (|2007) also performs well. 
In comparison, our algorithm produces somewhat more 
"bubbly" maps but appears to better capture the connec- 
tivity of HII regions. Both are an obvious improvement 
on directly filtering the linear density field. 

Of course, in our full algorithm we identify halos from 
the linear density field (rather than from simulations), so 
our method consumes comparable processing time to the 



one used to generate the leftmost panel in Figured! once 
the halos have been identified. Moreover, we are able 
to capture the "stochastic" component of the halo bias 
that causes the relatively large differences between the 
leftmost panel and the full RT simulation. That is, the 
algorithm used to generate the leftmost panel uses the 
large-scale linear dens ity fi el d to predict the distribution 
of halos (|Zahn et al.l I2005L l2007h . In reality, the rela- 
tion is not deterministic because of random fiuctuations 
in the small-scale modes comprising each region. This 
leads to nearly Poisson scatter in the halo number densi- 
ties (Sheth & Lemson 1999; Casas-Miranda et al. 2002|) 
that can substantially modify the bubble size distribution 
whenev er sources are rare, par ticularly early in reion- 
ization ([Furlanetto et al.l l2006al ). By directly sampling 
the small-scale modes to build the halo distribution, we 
better recover this scatter (at least statistically, as illus- 
trated by Fig. [2]). Another way to include this scatter 
is by directly sampling halos from an N-body simulation 
(as in IZahn et al.|[2QQ7l . or the second panel of Fig. [4j), 
although that obviously requires much more computing 
power. 

3.1. Ionization Maps 

Now that we have demonstrated in turn the success of 
our halo and bubble filtering procedures, we present the 
resulting ionization maps when the two are combined. 
In Figure El we show 100 Mpc x 100 Mpc x 0.5 Mpc 
slices through our 200^ ionization field at 2: = 10, 9, 
8.25, 7.25 {left to right across rows). With the assump- 
tion of C = 15.1, these redshifts correspond to xri =0.89, 
0.74, 0.53, 0.18, respecti vely. As has been pointed out by 
Furlanetto et ajj (|2004d ) , the neutral fraction is the more 
relevant descriptor; bubble morphologie s at a constant 
xhi vary little with redshift (see also iMcQuinn et al.l 



Fig. 4. — Slices from the ionization field at z = 6.89 created using different algorithms. All slices are 93.7 Mpc on a side and 0.37 Mpc 
deep, with the mean neutral fraction in the box being xhi = 0.49. Ionized regions are shown as white. The left-most panel was created 
by performing the bubble filtering procedure of Zahn et al. (2007) directly on the linear density field. The second panel was created by 
performing their bubble filtering procedure on their N-body halo field, but with the slightly different barrier definition in eq. (|14|) . The 
third p anel was created by performing our bubble filtering procedure described in §[3]on the same N-body halo field. The right-most panel 
ffrom Zahn et al.ll20Q7l ') was created using an approximate RT algorithm on the same halo field. 



l2006af ). The bottom- left panel corresponds to the halo 
field in the top-right panel of Fig. O generated on a 
high-resolution 1200^ grid. 

To quantify the ionization topology resulting from our 
method, we calculate the size distributions of both the 
ionized and neutral regions. We randomly choose a pixel 
of the desired phase (neutral or ionized), and record the 
distance from that pixel to a phase transition along a 
randomly chosen direction. We repeat this Monte Carlo 
procedure 10^ times. Volume- weighted probability dis- 
tribution functions (PDFs) produced thusly are shown 
by the solid curves in Figure [6] for ionized regions {top 
panel) and neutral regions {bottom panel). Curves corre- 
spond to (z, xm) = (10, 0.89), (9.25, 0.79), (8.50, 0.61), 
(8.00, 0.45), (7.50, 0.27), (7.00, 0.10), from left to right 
in the top panel, respectively (or from right to left in 
the bottom panel). All curves are normalized so that the 
probability density integrates to unity. 

It is useful to compare thes e distributions to the an - 
alytic bubble mass function of iFurlanetto et aH (|2QQ4c[ ): 
although this analytic approach is motivated by the same 
excursion set barriers as our semi-numerical approach, it 
does not account for the full geometry of sources. We 
compute the probability distribution from the analytic 
model by assuming purely spherical bubbles and convolv- 
ing with the volume-weighted distance to the sphere's 
edge: 



p{r) dr 



27rr^ dr 
(1 - Xm) 



J dRnb{R) 



1 



(15) 



where nb{R) is the comoving number density of bub- 
bles with radii betwee n R and R + dR (taken from 
IFurlanetto et al.ll20Q43 ). 

Several points are evident from Figures p and 
M As expected (e.g; IFurlanetto et all l2QQ4d l2QQ6al : 
iMcQuin n et al.l l2QQ6al ). there is a well-defined bubble 
scale at each neutral fraction, despite some scatter in 
the sizes. This scale also gets more pronounced (i.e. the 
PDF peaks more) as reionization progresses; this is a 
result of the cha ng:ing: shape of the und erlying matter 
power spectrum (Furlanetto et al."2006a'). 

Also, the purely analytic estimates underpredict the 
size distributions at all values of the neutral fraction, 
though they do become increasingly accurate as the neu- 
tral fraction decreases. This trend is perhaps counterin- 



tuitive, as the analytic model, which rests on the assump- 
tion of spherical bubbles, should perform best when the 
bubbles are isolated, as one would expect at earlier times, 
i.e. high neutral fractions. However, looking at the top- 
left panel of Fig. [5l the typical bubbles filling most of 
the ionized volume overlap due to the strong clustering 
of early sources and bubbles. This results in many "over- 
lapping pairs of spheres" at early times, resulting from 
merging HII bubbles sourced by clustered sources. Thus 
the spherical bubble-based analytic model underpredicts 
the true size distribution, using our "mean free path" def- 
inition of bubble si zes above. This e ffect was not noted 
by previous studies (|Zahn et al.ll2QQ7l ) , because they used 
a different definition of bubble sizes, based on spherical 
filters used to fiag regions in which xm < 0.1. As time 
progresses and the universe becomes more ionized, this 
"overlapping pair of spheres" effect becomes less and less 
dominant (see Fig. [5]), and the analytic model becomes 
increasingly more accurate. 

Finally, the size distributions of neutral regions pre- 
sented in the bottom panel of Fig. [6] are a new result and 
potentially important for the 21-cm signal (which origi- 
nates in neutral hydrogen, of course). In the later stages 
of reionization, when the topology has transformed to 
isolated neutral islands in a sea of ionized gas, this fig- 
ure pinpoints the typical sizes of "mostly neutral" pixels 
that continue to emit strongly. In contrast to the ionized 
regions, the neutral regions (defined in this way) do not 
grow substantially during reionization. From xm =0.89 
to Xm =0.1, the peak of the distribution shifts only by a 
factor of ~ 6, whereas the peak of the ionized region dis- 
tribution shifts by a factor of ~40 over the same range. 
The reason for this is also evident in Figure [5l even when 
the universe is mostly neutral, space is dotted with is- 
lands of ionized gas, such that our "mean free path" -type 
size distributions never become too large. The converse 
does not hold true for ionized regions. However, a slight 
parallel for ionized islands in a mostly neutral IGM, could 
be found in Lyman limit systems (LLS) inside larger HII 
regions fe.g. iBarka na fc Loebll2"QQ2l: iShapiro et al]l2QQ4l : 
Mirald a-Escude et al. 2000), though it is not clear how 
prevalent such neutral clumps are at high redshifts. 

Throughout this paper, we have used a L = 100 
Mpc "simulation" box. This size facilitates compar- 
ison of our results with those from recent hybrid N- 
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Fig. 5. — Slices through the 200"^ ionization field at z = 10, 9, 8.25, 7.25 {left to right across rows). With the assumption of = 15.1, 
these redshifts correspond to xui = 0.89, 0.74, 0.53, 0.18, respectively. All slices are 100 Mpc on a side and 0.5 Mpc deep. The bottom-left 
panel corresponds to the halo field in the top-right panel of Fig. [S] generated on a high-resolution 1200^ grid. 



body works (IZahn et all l2QQ7l: iMcQuinn et all l2QQ6al : 
llliev et al.l[2QQ6al : lTrac fc Cenll2QQ6r ). However, the speed 
of our semi-numerical approach allows us to explore 
larger cosmological scales while still consistently resolv- 
ing the small halos that could dominate the photon bud- 
get during reionization. As mentioned previously, exist- 
ing N-body codes must resort to merger-tree methods 
to populate their distributi on of small -mass halos, even 
for box sizes < 100 Mpc (iMcQuinn et al. 2006a). In 
this spirit, we present some preliminary results from a 
N = 1500^, L = 250 Mpc simulation, capable of directly 



resolving halos with masses M > 2.2 x lO^M©, with re- 
sulting mass functions accurate to better than a factor 
of two even at the smallest scale. This resolution pushes 
the RAM limit of our machine and so each redshift can 
take several hours to complete.^ 
In Figure [71 we compare size distributions of ionized 

^ We note here that our halo-finding algorithm requires signifi- 
cantly higher resolution than does predicting the ionization field 
directly from the linear density field smoothed on larger scales 
(Zahn et al. 2005, 2007). The latter method can be extended to 
even larger boxes, though at the price of a somewhat less accurate 
ionization map (compare the left and right panels in Fig. 0}. 
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Fig. 6. — Size distributions (see definition in text) of ionized 
(top panel) and neutral (bottom panel) regions. Curves correspond 
to (z, xui) = (10, 0.89), (9.25, 0.79), (8.50, 0.61), (8.00, 0.45), 
(7.50, 0.27), (7.00, 0.10), from left to right in the top panel, re- 
spectively (or from right to left in the bottom panel). Solid curves 
are produced from our simulation while dotted curves correspond 
to the analytic mass function. All curves are normalized so that 
the probability distribution integrates to unity. 



{top panel) and neutral {bottom panel) regions from our 
two different simulation boxes. Curves correspond to {z^ 
xm) = (9.00, 0.80), (8.00, 0.56), (7.00, 0.21), from left 
to right in the top panel, respectively (or from right to 
left in the bottom panel, respectively). Solid curves are 
generated from our fiducial, N = 1200^, L = 100 Mpc, 
simulation while dashed curves are generated from our 
larger simulation with N = 1500^, L = 250 Mpc. The 
cell size in all ionization maps is 0.5 Mpc on a side, with 
the efficiency parameter, adjusted to obtain matching 
values of xhi, and we set the minimum halo mass to 
^min = 2.2 X 10^ M0 even in the higher resolution runs 
for easier comparison. 

As reionization progresses, an increasing number of 
large HII regions are "missed" by the L = 100 Mpc sim- 
ulation. Interestingly, the analogous trend in the neutral 
region size distributions {bottom panel) is weaker. This is 
most likely because the "ionized island" effect limits the 
size distributions of neutral regions as described above. 

4. 21-CM TEMPERATURE FLUCTUATIONS 

A natural application of our "simulation" technique is 
to predict 21-cm brightness temperatures during reion- 
ization. The offset of the 21-cm brightness temperature 
from the CMB temperature, Xy, along a line of sight 
(LOS) at observed freq uency can be written as (e.g. 
iFurlanetto et al.ll2006bf ): 



sn{u)-- 



Ts 



Fig. 7. — Size distributions of ionized (top panel) and neutral 
(bottom panel) regions from different simulation boxes. Curves 
correspond to (z, xm) = (9.00, 0.80), (8.00, 0.56), (7.00, 0.21), 
from left to right in the top panel, respectively (or from right to 
left in the bottom panel). Solid curves are generated from our 
fiducial, = 1200^, L = 100 Mpc, simulation while dashed curves 
are generated from a larger simulation with A^ = 1500^, L = 250 
Mpc. The cell size in all ionization maps is 0.5 Mpc on a side, 
with the efficiency parameter, ^, adjusted to get matching values 
of Xm and the minimum halo mass set to Mmin = 2.2 x 10^ Mq 
for comparison purposes. 
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where Ts is the gas spin temperature, r^^ is the opti- 
cal depth at the 21-cm frequency vq^ ^ni is the physical 
overdensity (see discussion under eq. [Mj), H is the Hub- 
ble parameter, dVr/dr is the comoving gradient of the 
line of sight component of the comoving velocity, and all 
quantities are evaluated at redshift z = uq/u — 1. The 
final approximation makes the standard assumption that 
Ts ^ for al l redshifts of interest during reionization 
(e.g. lFurlanettoll2006h and also that dvr/dr <C H. We 
verify in our simulation that dvr/dr < H for all neutral 
pixels. 

Maps of (5T5(x, u) generated in this manner are shown 
in Figure [8l All slices are 100 Mpc on a side, 0.5 Mpc 
deep, and correspond to (z, xhi) = (9.00, 0.74), (8.25, 
0.53), (7.50, 0.27), from left to right. The top panels take 
into account the velocity correction term in eq. ([16]), 
while the bottom panels ignore it. 

As seen in Fig. [HI velocities typically increase the con- 
trast in temperature maps, making hot spots hotter and 
cool spots cooler. We also see that temperature hot 
spots, which correspond to dense pixels, tend to clus- 
ter around the edges of HII bubbles, especially smaller 
bubbles. This occurs because HII bubbles correlate with 
peaks of the density field and long- wavelength biases in 
the density field can extend beyond the edge of the ion- 
ized region. This enhanced contrast might be useful in 
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Fig. 8. — Brightness temperature of 21-cm radiation relative to the CMB temperature. All slices are 100 Mpc on a side, 0.5 Mpc deep, 
and correspond to (z, ^hi) = (9.00, 0.74), (8.25, 0.53), (7.50, 0.27), left to right. Top panels include the velocity correction term in eq. 
(|16|) , while the bottom panels do not. For animated versions of these pictures, see http: / / pantheon.yale.edu / ^am834 / Sim 



the detection of the boundaries of ionized regions with fu- 
ture 21-cm experiments. As reionization progresses most 
hot spots become swahowed up by HIT bubbles, and the 
effects of velocities diminish. 

In Figure [9] we plot the dimensionless 21- 
cm power spectrum, defined as A^iik^z) = 
ky{27r^V)_ (|fe(k,z)p)fe, where ^(x,^) = 
5Ti)((K^ z) / 5Ti){z) — 1. Solid blue curves take into 
account gas velocities, while dashed red curves do not. 
Curves correspond to (xhi, z) = (0.79, 9.25), (0.61, 
8.50), (0.45, 8.00), (0.27, 7.50), (0.10, 7.00), bottom to 
top. Error bars on the bottom dashed curve denote l-a 
Poisson uncertainties; fractional errors in a given bin 
are the same for all curves. As reionization progresses, 
small-scale power is traded for large-scale power, and the 
curves become fiatter. Note that, with our dimensionless 
definition of the power spectrum, curves with smaller 
xm have larger values of A2i(/c,^). This is because the 
mean brightness temperature offset drops quite rapidly 
as reionization progresses, since 6Ti){z) cx xri, but the 
scatter remains significant (see Fig. [6|) and thus the 
fractional perturbation, (52i(x, z), increases throughout 
reionization. 

Finally, in Figure [10] we plot dimensional power spec- 
tra, STb{z)^Al^{k,z). The curves correspond to (xhi, z) 
= (0.80, 9.00), (0.56, 8.00), (0.21, 7.00), top to bottom at 
large /c, respectively. The dotted green curves are gen- 
erated from a large, high-resolution "simulation", with 
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Fig. 9. — Dimensionless 21-cm power spectra for (xm, z) = 
(0.79, 9.25), (0.61, 8.50), (0.45, 8.00), (0.27, 7.50), (0.10, 7.00), 
bottom to top. Solid blue curves take into account gas velocities, 
while dashed red curves do not. 
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Fig. 10. — Dimensional 21-cm power spectra. The curves corre- 
spond to (xHi, z) = (0.80, 9.00), (0.56, 8.00), (0.21, 7.00), top to 
bottom at large k, respectively. The dotted green curves are gen- 
erated from a large, high-resolution "simulation", with = 1500"^ 
and L = 250 Mpc, with no velocity contribution to the power spec- 
tra. Solid blue curves and dashed red curves are generated with 
our fiducial = 1200"^ and L = 100 Mpc simulation, with and 
without the velocity contribution, respectively. 



N = 1500^ and L = 250 Mpc, with no velocity contribu- 
tion to the power spectra. Solid blue curves and dashed 
red curves are generated with our fiducial N = 1200^ 
and L = 100 Mpc simulation, with and without the ve- 
locity contribution, respectively. The cell size in all dT^ 
maps is 0.5 Mpc on a side, with the efficiency parame- 
ter, adjusted to achieve matching values of xm and 
the minimum halo mass set to Mmin = 2.2 x lO^M© for 
comparison purposes. 

As seen in Figures [9] and [TOl velocities make a mod- 
est contribution to the 21 cm power spectrum, boosting 
power on small scales early in reionization. Note that 
the apparent slight decrease in power at small k when 
velocities are included is well within the errors from av- 
eraging over the few modes available to us on the largest 
scales (e.g., see Poisson error bars on the bottom dashed 
curve in Fig. [9|). While the maximum ST^ value in our 
simulation box increases by a factor of a few when ve- 
locities are included, most of the pixels are only slightly 
affected. When the power spectrum is plotted in a di- 
mensional version, dT}){z)'^A2i{k, z), small-scale power is 
boosted by - 40% at (xri, z) = (0.80, 9.00), with this en- 
hancement monotonically decreasing as reionization pro- 
gresses. Linear theory predicts that velocities enhance 
the densi ty power spe ctrum by a factor of 1.87 when 
xhi = 1 (jKaiseriri987l ). In fact we do recover this en- 
hancement for a fully neutral IGM; however, as predicted 
by analytic models (McQuinn et al. 2006b), the ionized 
bubbles rapidly remove most of this amplification. 

Figure [To] also confirms the inferences drawn from Fig. 
[71 primarily that larger box sizes are needed to capture 
the ionization topology at the end stages of reionization. 
Comparing the dashed red to the dotted green curves in 
Fig. [TOl we note that our fiducial L = 100 Mpc simula- 
tions are accurate for scales smaller than A: > 0.2 Mpc~^ 
(or A < 30 Mpc). As reionization progresses, larger scales 
lose power more rapidly than in the L = 250 Mpc simu- 
lation. This is again evidence that very large scale simu- 
lations are needed to model the middle and late stages of 



reionization. Thus the speed and high resolution of our 
semi-numeric approach will be extremely useful in future 
modeling of reionization. 

5. CONCLUSIONS 

We introduce a method to construct semi-numeric sim- 
ulations that can efficiently generate realizations of halo 
distributions and ionization maps at high redshifts. Our 
procedure combines an excursion-set approach with ffist- 
order Lagrangian perturbation theory and operates di- 
rectly on the linear density and velocity fields. As such, 
our algorithm can exceed the dynamic range of exist- 
ing N-body codes by orders of magnitude. As this is 
the main limiting factor in simulating the ionized bubble 
topology throughout reionization, when ionized regions 
reach scales of tens of comoving Mpc, this will be partic- 
ularly useful in such studies. Moreover, the efficiency of 
the algorithm will allow us to explore the large parame- 
ter space required by the many uncertainties associated 
with high-redshift galaxy formation. 

We find that our halo finding algorithm compares well 
with N-body simulations on the statistical level, yield- 
ing both accurate mass functions and power spectra. We 
have not yet compared our halo distribution with sim- 
ulations on a point-by-point basis, but we do not ex- 
pect perfect agreement because of the vagaries of the ex- 
cursion set approach. However, it is encouraging that 
a very similar algorithm independently developed by 
iBond fc MversI (|l996al ) fares quite well in a comparison 
of high-mass halos. 

Our HIT bubble finding algorithm captures the bubble 
topology quite well, as compared to ionization maps from 
ray-tracing RT algorithms at an identical xm- Our al- 
gorithm is similar to other codes, although we build the 
ionization map from our excursion set halo field rather 
than directly from the linear de nsity field o r fro m halo s 
found in an N-body simulation (|Zahn et al.ll2005L I2007D . 
Compared to codes built only from the linear density 
field, we can better track the "stochastic" component of 
the bias, though at the cost of somewhat more compu- 
tation and a harder limit on resolution. On the other 
hand, our scheme is much faster than using an N-body 
code and offers superior dynamic range. 

We create ionization maps using a simple efficiency pa- 
rameter and compute the size distributions of ionized 
and neutral regions. Our size distributions are gener- 
ally shifted to larger scales when compare d with purely 
analytic models (Furlanetto et al.l l2004cf ) at the same 
mean neutral fraction. The discrepancy lies in the fact 
that, at their core, the purely analytic models are based 
on ensemble- averaged distributions of isolated spheres. 
Hence they do not capture overlapping bubble shapes, 
which are most important at large xm (when the bubbles 
are small and random fluctuations in the source densities, 
as well as clustering, are most important). 

In this paper, we have confined ourselves to a sim- 
ple ion ization c riterion (essentially photon counting; 
iFurlanetto et aP l2004cl ). However, our algorithm can 
easily accommodate more sophisticated prescriptions, so 
long as they ca n be expressed either w ith the excursion 
set for malism (|Fur lanetto & Oh 20Q5|; IFurlanetto et al.l 
'2006a') or built from the halo field (in a similar way to 
semi-analytic models of galaxy formation embedded in 
numerical simulations). 
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We also use our procedure to generate maps and power 
spectra of the 21-cm brightness temperature fluctuations 
during reionization. We note that temperature hot spots 
generally cluster around HII bubbles, especially in the 
early phases of reionization. Because HII bubbles cor- 
relate with peaks of the density field, long- wavelength 
biases in the density field can extend beyond the edge of 
the ionized region, with the resulting overdensities ap- 
pearing as hot spots. This effect might be useful for 
detecting the boundaries of ionized regions with future 
21-cm experiments. We study the imprint of gas bulk ve- 
locities on 21-cm maps and power spectra, an effect which 
was not included in previous studies. We find that ve- 
locities do not have a major impact during reionization, 
although they do increase the contrast in temperature 
maps, making some hot spots hotter and some cool spots 
cooler. Velocities also increase small-scale power, though 
the effect decreases with decreasing xm- 

We also include some preliminary results from a sim- 
ulation run with the largest dynamical range to date: 
a 250 Mpc box which resolves halos with masses M > 
2.2 X 10^ Mq. This simulation run confirms that ex- 



tremely large scales are required to model the late stages 
of reionization, xui ^0-5, when the typical scale of ion- 
ized bubbles becomes several tens of Mpc. 

The speed and dynamic range provided by our semi- 
numeric approach will be extremely useful in the mod- 
eling of early structure formation and reionization. Our 
ionization maps can be efficiently folded into analyses of 
current and upcoming high-redshift observations, espe- 
cially 21-cm surveys. 
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